#include <stdio.h>
#include <stdlib.h>

#include <math.h>
#include <mpi.h>

#include "fft_parallel.h"

#define SIZE 524288
#define PI (M_PI*2.0)

#define START_BCAST_EVENT 1
#define END_BCAST_EVENT 2
#define START_SEND_EVENT 3
#define END_SEND_EVENT 4
#define START_RECV_EVENT 5
#define END_RECV_EVENT 6

int main(int ac, char **av)
{
    FILE *fin;
    int ix, nproc, myrank, myid;
    int count = 1;
    int length = 30;
    char name[length];

    double *data = (double *) calloc(SIZE, sizeof(double));
    double *data_temp = (double *) calloc(SIZE, sizeof(double));
    double *data_orig = (double *) calloc(SIZE, sizeof(double));

    MPI_Init(&ac, &av);
    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Get_processor_name(name, &length);
    MPI_Status status;

    for (ix = 0; ix < SIZE; ++ix) {
        data[ix] = 0.0;
        data_temp[ix] = 0.0;
    }

    if(myrank == 0) {
    	fprintf(stderr, "Odczyt danych... ");
        fin = fopen("dane.txt", "r");
	ix = 0;
        while(fscanf(fin, "%lf", &data_orig[ix++]) != EOF) ;

	fclose(fin);
	fprintf(stderr, "ok\n");

        for (ix = 0; ix < SIZE; ++ix)
            data[ix] = data_orig[ix];

	fprintf(stderr, "Algorytm BRP... ");
        brp(data);
	fprintf(stderr, "ok\n");

	MPI_Bcast(data, SIZE, MPI_DOUBLE, 0, MPI_COMM_WORLD);

        fft(data, nproc, myrank);

        while (count < nproc) {
            MPI_Recv(data_temp, SIZE, MPI_DOUBLE, MPI_ANY_SOURCE, 2, MPI_COMM_WORLD, &status);
            for (ix = 0; ix < SIZE; ++ix) {
                    if(data_temp[ix] != 0)
                        data[ix] = data_temp[ix];
            }
            ++count;
        }
    } else {
        MPI_Bcast(data, SIZE, MPI_DOUBLE, 0, MPI_COMM_WORLD);

        fft(data, nproc, myrank);

        MPI_Send(data, SIZE, MPI_DOUBLE, 0, 2, MPI_COMM_WORLD);

    }

    if (myrank == 0) {
        for (ix = 0; ix < SIZE; ++ix) {
            printf("%f\n", data[ix]);
        }
    }
    MPI_Finalize();

    return 0;
}

void brp(double *data)
{
    int brp_index[SIZE];
    int i, m, x;
    double tmp[SIZE];

    data[0] = 0;
    x = 0;
    for(i = 1; i < SIZE; ++i) {
        m = SIZE;
        do {
            m = m >> 1;
            x = x ^ m;
        } while (x < m);
        brp_index[i] = x;
    }

    for(i = 0; i < SIZE; ++i)
	tmp[i] = data[i];

    for(i = 0; i < SIZE; ++i)
	data[i] = tmp[brp_index[i]];
}

void fft(double *data, int nproc, int myrank)
{
    int d, i, j, k, s;
    double data_img[SIZE];
    double tr, ti, ur, ui;
    double wmr, wmi, wr, wi;

    for (i = 0; i < SIZE; ++i)
        data_img[i] = 0;

    for (s = 0; s < 16; ++s) {
        d = pow(2, s);
        wmr = cos(PI / (double)d);
        wmi = sin(PI / (double)d);
        for (k = 0; k < SIZE; k += d) {
            wr = 1;
            wi = 0;
            for (j=0; j<(d/2); j++) {
                if ( (j + k)%nproc == myrank || (j + k + (d / 2))%nproc == myrank )
                {
                    tr = wr*data[k + j + (d/2)] - wi*data_img[k + j + (d/2)];
                    ti = wr*data_img[k + j + (d/2)] + wi*data[k + j + (d/2)];
                    ur = data[k + j];
                    ui = data_img[k + j];
                    data[k + j] = ur + tr;
                    data_img[k + j] = ui + ti;
                    data[k + j + (d/2)] = ur - tr;
                    data_img[k + j + (d/2)] = ui - ti;
                }
                wr = wr * wmr - wi * wmi;
                wi = wr * wmi + wi * wmr;
            }
        }
    }

    for (i = 0; i < SIZE; ++i)
        data[i] = sqrt( pow(data[i], 2) + pow(data_img[i], 2) );
}
